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^ . Abstract 

*^ ■ This paper proposes a concise coding of the cehs of Ji-dimensional finite regular grids. 

It induces a simple, generic and efficient framework for implementing classical dig- 
ital topology data structures and algorithms. Discrete subsets of multidimensional 
images (e.g. regions, digital surfaces, cubical cell complexes) have then a common 
and compact representation. Moreover, algorithms have a straightforward and effi- 
tyj \ cient implementation, which is independent from the dimension or sizes of digital 

images. We illustrate that point with generic hypersurface boundary extraction al- 
gorithms by scanning or tracking. This framework has been implemented and basic 
operations as well as the presented applications have been benchmarked. 
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1 Introduction 



Many applications in the image analysis field need to represent and manipulate 
discrete subsets of digital spaces. As the image data become larger, the data 
^ ■ structures required to represent these sets should be as compact as possible. 

■ Moreover, algorithms designed on these structures should be not only efficient 

theoretically, but also efficient in practice. Algorithms defined formally should 
have a straightforward implementation. Last but not least, 3D and 4D image 
datasets are now more and more common. It becomes necessary to have a 
unified framework for programming applications dealing with n-dimensional 
data. By this way, algorithms are both generically defined and implemented. 

There exist several approaches to define the topology of multidimensional reg- 
ular digital spaces: (1) adjacency graphs as pioneered by Rosenfeld, (2) ori- 
ented graphs as proposed by Herman [S], (3) cellular complexes as proposed 
by Kovalevsky [11], or equivalently Khalimsky's spaces [7] and interpixel rep- 
resentations. This paper deals mostly with the third approach, although our 
framework can express either approaches (the two first approaches manipulate 
a restricted set of the elements defined in the spaces of the third approach) . 
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The cellular decomposition of the Euclidean n-dimensional space into a 
regular grid forms a cellular complex C". This structure has been introduced 
in digital topology by Kovalevsky [TT] for 2D and 3D applications. It has been 
shown that the topology induced on C" is equivalent to a digital topology K", 
sometimes called Khalimsky topology [71I8] . Many authors have explored the 
theoretical properties of the space C" (or equivalently K") [7|l8|[TT] applications 
[2][3]|9][T0l[T^ . These works show that the definition of consistent high level data 
structure over images relies on a low-level representation which is the regular 
cellular decomposition of the image support. It is thus critical to represent 
efficiently arbitrary cells of C", small subsets of C" and specific subsets of 
C" (e.g. complexes, digital surfaces). However, the litterature does not refiect 
this observation. Indeed, spels are often coded with an array of coordinates, 
surfels as pairs of adjacent spels or a spel with a direction, other cells are 
generally implicitly represented. Consequently, storing elements or subsets of 
C" is cumbersome; algorithms are frequently rewritten at the implementation 
stage to avoid any reference to non elementary kinds of cells. 

In this paper, we choose another approach, which is first to show how to repre- 
sent an arbitrary cell of C" with a binary cell code and secondly to design data 
structures over this representation. Because of the regularity of C", each cell 
code holds all the information on the cell: the cell topology (dimension, open 
or closed along a coordinate, adjacent and incident cells) and geometry (co- 
ordinates in Z", centroid, normal and tangent vectors) can be computed from 
the code without any other information. The proposed framework is suited 
both to formal representation and proofs and to straightforward implementa- 
tion in a programming language. The paper is organized as follows: (i) coding 
of cells and implementation of low-level digital topology definitions, (ii) defi- 
nition of data structures for subsets of (e.g., digital surfaces, complexes), 
(iii) application to digital boundaries extraction in multidimensional images. 
We emphasize that all operations and algorithms have the same definitions 
and implementation whatever the dimension of the space. All experiments 
and benchmarks presented in this paper were made on a PC with a Celeron 
SOOMhz processor, 128Mb of memory, 128Kb of cache (which is a basic work- 
station). The proposed framework was implemented in C-I--I-. Due to limited 
space, the reader is referred to [12] for more details. 



2 Coding cells of digital spaces C" 

2.1 Cellular decomposition C"; binary coding of unoriented cells 

We denote by C" the set of parts of the n- dimensional Euclidean space M" 
such that c G C" is equivalent to c = /i x . . . x where Jj is a subset of M of 
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the form {zi} or ]zi, Zi + 1[ with Zi G Z. The complex C" is a partition of M". 
We call k-cell an element c G C" such that c has A; /j of the form ]zi, Zi + 1[ 
(and therefore (n — k) li of the form {^j}). The dimension of c is A;. The closure 
Cl(c) of a cell c is the set of cells c' of C" which have the following form: (i) on 
coordinate where c is a point {zi}, the cell c' must also be the same point, (ii) 
on coordinate where c is an open segment + 1[, c' can be either the same 
open segment or the point {zi} or the point {zi + 1}. The open star Op(c) 
of a cell c is the set of cells c' of such that c' G Op(c) -v^ c G Cl(c'). The 
bounding relation < between two cells c and c' is then defined as c < c' iff 
c G Cl(c') \ c'. With these definitions, the set C" equipped with the dimension 
mapping and the bounding relation is a cellular complex. A cubical cell complex 
K is then defined as any set of cells in a finite image. The dimension of K is 
the maximum of the dimensions of its cells. Open stars and closure of cells in 
a complex K are defined naturally. 



It is clear that the elements of C" represent low-level elements of n-dimensional 
digital images: the spels (pixels in 2D and voxels in 3D) are the n-cells, the 
(unoriented) surfels (a pair of adjacent spels) are the ri — 1-cells, the vertices of 
the spels and of the surfels, or pointels, are the 0-cells. An object is then a set of 
?7,-cells, a digital surface is a set of n — 1-cells (oriented or not, see Section [2^2]) . 
a curve is a set of connected 1-cells. Therefore, all classical subsets of digital 
spaces have a natural definition as specific subsets of C". From now on, we 
will assume that we are working in a finite n-dimensional image forming a 
parallelepiped in Z". We denote by M* the inclusive upper bound for the ^-th 
coordinate of any spel. All coordinates have as lower bound. 



As shown by Kong et al. [8] , the topology of C" is equivalent to the topology 
of the Khalimsky digital space K", which is the cartesian product of n con- 
nected ordered topological spaces (COTS). A COTS can be seen as a set of 
ordered discrete points, like Z, whose topology alternates closed points and 
open points. If we define even points of Z as closed and odd points of Z as 
open, each point of K" is then identified by its n integer coordinates, whose 
parities define its topological properties. 



Consequently, any Ac-cell c of C" has exactly one corresponding point in 
with coordinates {x%;, . . . ,x^~^). We propose to code c as one binary word 
code(c) = a . . . . . . x° , called the unsigned code of c, as follows: 



• The i-th coordinate is coded by its binary decomposition after a right- 
shift ( div 2). We say that is the i-th digital coordinate of c. 

• All coordinates are packed as one binary word (from x"~^ to x°). Every 
coordinate is allocated a fixed number of bits A^, given by Aj = log2(M*) -|- 1. 

• The parity of all coordinates are also packed as an n-bits word a with 
a = J^ii-C^K ^od 2)2*. a is called the topology of c. 
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Spels have a topology word composed of I's, whereas pointels have a topology 
word made of O's. Surfels have one and n—1 I's in their topology word. The 
coordinate where a surfel c has a in its topology word is called the coordinate 
orthogonal to the surfel c and is denoted by ± (c). This coding implies that any 
cell of finite digital images can be coded as an integer number with fixed size. 
Any register of a processor may thus store a cell if the image is not too bigQ 
We define the adjacency between cells independently of the cell topology. 

Definition 1 Two cells p andq with topo(p) = topo(g) are /-adjacent if their 
respective coordinates jo* and (f differ by at most 1 and if the infinite norm of 
the vector — q"'~^, . . . — g°) is no more than I. 

The 1-adjacency thus defines the 4-adjacency (resp. 6-adjacency) on pixels in 
2D (resp. on voxels in 3D) and the 2-adjacency defines the 8-adjacency (resp. 
18-adjacency) on pixels in 2D (resp. on voxels in 3D). We define the incidence 
relation as below. The proposition that follows shows that all the topological 
structure of C" can be obtained from the incidence relation. 
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. A cell p is low incident (resp. up incident j to a cell 
q if there is a sequence of cells Cq = p, Ci, . . . , Ck = q such that Cj is low 
1-incident (resp. up 1-incident) to Cj+i. 

Proposition 3 The set of cells low incident to a cell c is equal to Cl(c) \ c. 
The set of cells up incident to a is equal to Op(c) \ c. 

Figure [H summarizes the number of elementary operations necessary to execute 
basic cell operations. Their implementation is fully generic. We have bench- 
marked these operations and the results show that cell codes compete with 
statically defined structures (e.g. fixed size arrays) and are much faster than 
dynamically allocated structures, classically used for generic programming. 



2.2 Oriented cells, boundary operators, bels, boundary of an object 



In some applications, it is convenient to orient the cells (as positive or neg- 
ative). For instance, digital surfaces as proposed by Herman and Udupa are 

^ 32 bits are sufficient to code every cell of a 1024 x 1024 x 512 3D image, which is 
more than enough for current biomedical applications. 
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Fig. 1. Number of elementary operations needed to perform the fohowing tasks: (i) 
coding a vector of n Khahmsky coordinates as a cell, (ii) getting the topology or 
one coordinate of a cell, (iii) comparing if two cells are identical, (iv) setting the 
coordinate of a cell, (v) computing a 1-adjacent cell, (vi) checking if two cells are 
/-adjacent, (vii) computing a 1-incident cell, (viii) checking if two cells are /-incident. 

composed of oriented pairs of voxels: one voxel is in the interior of the surface, 
the other in the exterior. Orienting a surfel means in this case to define where 
are the interior and exterior voxels 1-up-incident to the surfel. Digital surface 
tracking algorithms rely on this orientation for a consistent output. Classi- 
cal combinatorial topology associates an orientation to each cell of a cellular 
complex. Oriented cells are then useful to implement boundary operators over 
complexes and to compute topological invariants. 

We therefore define the signed code of a cell c with orientation bit s (0 is 
positive, 1 is negative) by adding the bit s between the topology a of c and its 
digital coordinates follows: 
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s 













. The opposite cell — c 

of c is the same cell as c but with opposite sign. Boundary operators, which 
can be seen informally as an oriented version of incidence, are essential in 
combinatorial topology : for instance, they define the topology of polyhedral 
complexes. We have now to "orient" the incidence relation. 
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, is called the lower boundary of the 
c along coordinate ij. The lower boundary Ac oj c is then the set of cells 



U/=o,...,fcAj,c. 



The lower boundary of c thus corresponds to the set of cells 1-low-incident to c 
with specific orientations. The upper boundary V of a cell is defined symmetri- 
cally (the upper boundary is taken on topology bits set to 0). It can be shown 
that this definition of boundary operators induces that any cubical cell com- 
plex is a polyhedral complex. Although this is outside the scope of this paper. 
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boundary operators along with chains are used to define the (co)homology 
groups of complexes, which are well known topological invariants. This topo- 
logical tool is readily applicable to cubical cell complexes. 

In the remainder of the paper, the set O is an object of the image / with 
an empty intersection with the border of /. Assume that all spels of O are 
oriented positively. We merge the sets Ap with p & O with the rule that two 
identical cells except for their orientation cancel each other. The resulting set 
of oriented surfels is called the boundary of O, denoted by dO. It is a digital 
surface, whose elements are called bels of O. The following result states that 
the boundary of O is indeed the digital surface separating the spels of O from 
the spels of the complement of O. 

Proposition 5 Let c be a bel of dO. Then Vc contains two spels: one posi- 
tively oriented and belonging to O, the other negatively oriented and not be- 
longing to O. Any path of 1-adjacent spels from an element of O to an element 
not in O crosses dO. 

One way to compute the digital surface bordering a set of spels O is by apply- 
ing the lower boundary operator on each element of O and removing oppositely 
oriented identical cells. The time complexity of this algorithm is thus linear 
with the number of spels of O. 

2.3 Followers of surf el, bel adjacency, digital surface tracking 

The bel adjacency defines the connectedness relations between bels bounding 
an object. It has two nice consequences: (i) the boundary of an object can be 
extracted by tracking the bels throughout their bel adjacencies [iP] : (ii) sets of 
surfels can be considered as classical Euclidean surfaces, where one can move 
on the surface in different orthogonal directions (2 in 3D). The second reason 
is essential for defining the geometry of digital surfaces [12]. We present here 
a definition of bel adjacencies that is essentially equivalent to the definition 
of [6], but easier to implement in our framework. We start by defining which 
surfels are potentially adjacent to a given bel with the notion of follower. We 
then define two kinds of bel adjacency for each pair of coordinates. 

Definition 6 We say that an r-cell q is a direct follower of an r-cell p, p ^ 
±q, if Ap and Aq have a common r — 1-cell, called the direct link from p to q, 
such that this cell is positively oriented in Ap and negatively oriented in Aq. 
The cell p is then an indirect follower of q. 

It is easy to check that any surfel has 3 direct followers and 3 indirect followers 
along all coordinates except the one orthogonal to the surfel. We order the 
followers consistently for digital surface tracking (see Figure [2]). 
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eA,-q (3) 
e V, + 6' (2) 
eA^+p (1) 

Fig. 2. Direct followers of a surfel b along coordinate x. 



Definition 7 Let b be an oriented n — 1-cell, such that V& = {+p, —q}- Let 
j be a coordinate with j 7^_L (p). The three direct followers of p along j are 
ordered as follows: (1) the first direct follower belongs to /S.j+p, (2) the second 
direct follower belongs to V j + b' with -\-b' direct link in Ajb, (3) the third direct 
follower belongs to Aj — q. 



Intuitively, when tracking a digital surface, you have 3 different possibilities 
for a move along a given coordinate. This is true for arbitrary dimension. The 
following definition shows which one to choose at each step . It is in agreement 
with the definitions of bel adjacencies proposed by Udupa [T5] . 

Definition 8 Let b be a bel of dO, such that V6 = {+p, —q} (thus p e O 
and q ^ O). For any coordinate j 7^_L [b), the bel b has one interior direct 
adjacent bel (resp. exterior direct adjacent bel^ which is the first (resp. last) 
of the three ordered direct followers of b along coordinate j that is a bel of O. 
The bel adjacency is the symmetric closure of the direct bel adjacency. 



In 3D, the interior (resp. exterior) bel adjacency along all coordinates induces 
the classical (6,18) bel-adjacency (resp. (18,6) bel-adjacency). Interior and 
exterior bel adjacencies can be mixed for different coordinate pairs. This might 
be useful in an application where the image data are not isotropic (e.g., some 
CT scan images, confocal microscopy). Computing the bel adjacent to a given 
one is very fast since it required [12]: < 11 binary or integer operations, < 3 
shifts, < 14 lut accesses, < 9 conditional tests, and 1 or 2 "is in set" operations. 
The next section will show that the "is in set" operation can be done in four 
elementary operations. 



The following theorem, which comes from the fact that cubical cell complexes 
are polyhedral complexes, is interesting to speed up digital surface tracking 
algorithm: as its corollary, tracking only direct adjacent bels is sufficient to 
extract the whole digital surface component that contains the seed bel. It is 
more complex to show that bel components correspond to interior and exterior 
components of spels (see [T3IT5] where this is proven for some bel adjacency 
relations) . 

Theorem 9 The transitive closure of the direct bel adjacency from a bel b of 
do coincides with the transitive closure of the bel adjacency from b. 
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Table 1 

This table shows some properties of classical set data structures. The symbol + 
indicates that it is only amortized time complexity. 

set structure STL class is in set ? other set ops memory (bytes) 

dynamic array vector 0{m) 0{m)+ ~ 4m 

linked list list 0{m) 0{m) ~ 24m 

RB-tree set 0(logm)+ 0(logm)+ ~ 32m 

hashtable hash_set 0{l)+ 0(1)+ w 4?n' + 20m 

3 Data structures built over cells 



Since any kind of cell is coded as an integer number, data structures coding 
sets of cells are easily derived from standard data types. Table [1] displays 
the traditional set data structures, their implementation in C++ as ternplate 
classes, the time complexities of some operations, and the memory cost] ^ I 

These data structures are adapted to sets of cells of reasonnable size. Very 
small sets should be defined as vectors, lists may be used to represent 
medium size contours. Other medium size sets should be represented with 
sets or hash_sets. If the hash_set seems rather efficient for most operations 
(at least from an asymptotic point of view), it is memory costly: 28 Mbytes are 
necessary to represent a digital surface with 1,000,000 bels (and m' = 2m). 
Moreover the memory is very fragmented and the cache is thus not efficient. 
As it is shown later on digital hypersurface tracking algorithms, amortized 
constant time does not mean very fast. 

We present another data structure to represent a set of cells, which exploits 
the properties of the cell coding. The size of the data structure is dependent 
only on the size of the image. The time complexity of all operations is then 
independent from the number of cells represented. This data structure, called 
the CharSet, is a characteristic function that assigns one bit to each cell of 
the space. Since we will often manipulate sets of cells that contains specific 
kinds of cells (e.g.. a digital surface is made of surfels), we present two ways 
to define this structure. 

Definition 10 A MinCharSet is an array tbl of s hits, where s is one plus the 
difference between the highest possible cell code MAX and the smallest possible 
cell code min. Selecting the bit characteristic of the presence of a given cell c 
is done with tbl [(c-min)»5]&(l«(c&0xlf )) for 32-bits words. 

Definition 11 ^ LUTCharSet is an array tbl of s bits and a look-up ta- 
ble lut, where s and lut are dependent on the set of cells (see TableWi). 

^ A cell is stored in a 32-bits word. We suppose that 12 bytes are necessary to store 
information about one dynamically allocated memory area (e.g. holds for Linux) 
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Table 2 

This table defines the way LUTCharSets store various specific sets of cells. 

size s 256^ image 
(bits) size (Mb) 



set of cells topologies a lut(a) 



set of 
spels 



set of 

oriented 

r-cells 



1...1 





Oil 


. . 1 


unoriented 


101 


. . 1 


digital 






surface 


1... 


110 



0...0011...1 



0...0101...1 



1...1100...0 



0...0 



E bits 






0...0 


1 


0. . .0 



re — 1 



0. . .0 



bits 



1 



0...0 



re. 






0...0 


1 


0...0 



2(;!)2E^» 12(r = l,2) 



1 + bits 



Selecting the bit characteristic of the presence of a given cell c is done 
with tbl [ (lut [tope (c) ] +sign_coords (c) ) »5] & (1« (c&Oxlf ) ) for 32-hits 
words. 

The LUTCharSet is more compact than the MINCharSet for some sets of cells 
(and the higher the dimension the more it is) but the access to the character- 
istic bit of a cell is a bit slower It is now clear why the bit defining the sign 
of an oriented cell is inserted between the topology and the coordinates of the 
cell: with this coding, both CharSets use exactly twice more memory for sets 
of signed cells compared with sets of unsigned cells. 

Knowing if a cell belongs to a CharSet or any other atomic set operation 
(add/remove an element) are 0(1) operations. All global set operations (like 
union, intersection, difference, complement) between CharSets are imple- 
mented as standard bit operations between arrays of binary words. Their time 
complexities are linear in the size of the array. Moreover, the implementation 
of set operations for any set of cells (arbitrary dimension, set of spels, oriented 
digital surface, set of r-cells, etc) is done only once as bit operations between 
arrays of binary words. To give an idea of the efficiency of this representation, 
inverting a set of spels defined in a 512^ image takes 0.40s (134, 217, 728 spels, 
3ns per spel), difference between two sets of spels in the same image takes 
0.80s. Furthermore, an unoriented digital surface in a 256^ image can hold up 
to 50,331,648 surfels for a 8Mb memory cost (or 6Mb for LUTCharSet). To 
conclude this section, unsigned sets are twice less costly to store. They should 
be used when possible. For instance, digital surfaces that are boundaries of a 
set of spels are always orientable surface. Digital surface tracking can thus be 
done with unoriented digital surfaces. 
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// Track (B) algorithm. 
// do must he dosed. 
CharSet 

Space: :track( CharSet O, Cell b, 

BelAdj A ) 
{CharSet S = emptySurfelSet(); 
queue<Cell> L; // queue of bels 
L.push( b ); // starting bel 
while ( ! L.emptyO ) { 
Cell p = L.popO; // current bel 
/ / On all coord where p open 
for ( int j = 0; j < dim(); ++j ) 
if (j != orthDir(p)) { 
/ / Track direct followers 
Cell q = A.directAdj(0,p,j); 
if (! S.islnSet(q)) { 
S.add( q ); 
L.push( q ); 

} 

} 

} 

return S; 

} 



// Track (C) algorithm. 
// do must be closed. 
CharSet 

Space:: track ( CharSet O, Cell b, 

BelAdj A ) 
{CharSet S = emptySurfelSet(); 
queue<Cell> L; // queue of bels 
hst<Cell> T; // "tail" of bdry 
L.push( b ); // starting bel 
T.multiplelnsert( b, dim() - 1 ); 
while ( ! L.emptyO ) { 
Cell p = L.popO; // current bel 
for ( int j = 0; j < dim(); ++j ) 
if ( j != orthDir( p ) ) { 
Cell q = A.directAdj( O, p, j ); 
if ( T.find( q ) ) // already 

T.remove( q ); // extracted 
else { 
S.add( q ); L.push( q ); 
T. mult iplelnsert ( q, dim ( ) - 2 ) ; 

} }} 

return S; // T is empty at loop end 

} 



Fig. 3. Two digital hypersurface tracking algorithm: the Track (B) algorithm re- 
quires an efficient "is in set" operation, the Track (C) algorithm stores the hst of 
cehs that will be hit again by the tracking. For the set T in (b), we have tried both 
list and multiset. The former was much faster than the later in our experiments. 



4 Digital boundary extraction by scanning and tracking 



We have implemented several digital hyper-surface extraction algorithms which 
build the digital surface that is the boundary of a given object O. Scanning 
algorithms examine every spel neighborhood to detect the presence of a bel. 
They only require the set O as input. Digital surface tracking algorithms re- 
quire an initial bel b and a bel adjacency A to extract the component of the 
boundary of O that contains b. As described in Section 12.31 defining the bel 
adjacency A is deciding for each pair of coordinates whether A is interior or 
exterior along this plane. Figure [3] shows how to write generic digital surface 
tracking algorithms with our framework. The implementation in C++ is very 
close to the formal specification of the algorithm (see [B]). 

In the experiments, the object O was a digital volumic ball. Table [3] lists the 
running times necessary to extract 80 for balls of various radii and dimen- 
sions. The Scan (A) algorithm scans the whole image to find boundaries. The 
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Table 3 



Running times for several boundary extraction algorithms (see text). 



Space 
size 


xtad. 


Nb spels 


Nb surf. 


Scan 
(A) 


Scan 
(B) 


Track 
(A) 


Track 
(B) 


Track 
(C) 


40962 


2000 


12566345 


16004 


2.07s 


2.00s 


< 0.01s 


< 0.01s 


0.01s 


128^ 


30 


113081 


16926 


0.38s 


0.03s 


0.01s 


0.01s 


0.06s 


128^ 


60 


904089 


67734 


0.39s 


0.34s 


0.07s 


0.05s 


0.57s 


256^ 


120 


7236577 


271350 


3.15s 


2.70s 


0.36s 


0.32s 


5.24s 


5123 


240 


57902533 


1085502 


25.1s 


21.2s 


1.88s 


1.85s 


50.6s 


64^ 


30 


4000425 


904648 


4.26s 


4.00s 


1.91s 


1.37s 


4748s 



Scan (B) algorithm scans the parallelepipedic subspace containing the balL 
The Track (A) algorithm extracts open or closed boundaries from a starting 
bel (it follows both direct and indirect bel adjacencies). Track (B) and (C) 
algorithms extract only closed boundaries from a starting bel (they follow 
only direct bel adjacencies). All these algorithms are written generically and 
make no assumption on the dimension of the image. The benchmarks show 
that scanning algorithms depend on the size of the scanned subspace and that 
tracking algorithms depend on the number of surfels in dO. Running times 
are excellent since each bel is tracked in k, 1.7/xs in 3D (and ~ 1.5/is in 4D). 
Note that Track (B) algorithm is much faster than Track (C) algorithm. This 
is because CharSets are efficient for the query "is a cell in a given set ?" . 



5 Conclusion 

We have presented a binary coding of every cell of the digital space C". This 
coding contains all the topological and geometric information on the cell. It 
allows the design and implementation of generic low- level algorithms that deals 
with subsets of C". Compact and efficient data structures can be built with this 
coding. We illustrated the potential of this framework with a classical digital 
topology application: boundary extraction. Arbitrary dimensional algorithms 
are readily implemented in this framework and benchmarks have proved that 
the resulting code is surprisingly efficient in practice. Other digital topology 
and geometry applications may be found in |12j . 
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